hw3_1

library(lattice)
library(MASS)
library(e1071)
library(latticeExtra)

## Loading required package: RColorBrewer

library(corrplot)

concr.data <- read.csv("data/Concrete_Data.csv", comment.char = "#")
panel <- function(...) {  
  panel.xyplot(...)
  panel.loess(...)
}

xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$Superplasticizer, panel=panel)

xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$Cement, panel=panel)

Видна отличная линейная зависимость прочности от цемента, а также от Superplasticizer(далее Суперпластик). Ну это и понятно, т.к. они являются главными скрепляющими и упрочняющими ингредиентами

xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$Water, panel=panel)

С водой происходит не очень понятный скачек вверх после систематического опускания. Вообще говоря мы достигаем минимума при самой большой плотности точек, а далее точек совсем мало, и возможно, дальшейшие результаты можно сбросить на недостаточное количество данных. Попробуем его позже разобрать подробнее.

xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$CoarseAggregate, panel=panel)

xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$FineAggregate, panel=panel)

Аггрегаты очень похожы, и переваливая за 810 практически не влияют на прочность, скорее всего они будут малозначимыми

xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$Age, panel=panel)

xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$BlastFurnaceSlag, panel=panel)

xyplot(concr.data$ConcreteCompressiveStrength ~ concr.data$FlyAsh, panel=panel)

Возраст выглядит совсем неоднородным, видимо его нужно будет сделать фактором

model1 <- lm(ConcreteCompressiveStrength ~ (.), data=concr.data)
summary(model1)

## 
## Call:
## lm(formula = ConcreteCompressiveStrength ~ (.), data = concr.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.654  -6.302   0.703   6.569  34.450 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -23.331214  26.585504  -0.878 0.380372    
## Cement             0.119804   0.008489  14.113  < 2e-16 ***
## BlastFurnaceSlag   0.103866   0.010136  10.247  < 2e-16 ***
## FlyAsh             0.087934   0.012583   6.988 5.02e-12 ***
## Water             -0.149918   0.040177  -3.731 0.000201 ***
## Superplasticizer   0.292225   0.093424   3.128 0.001810 ** 
## CoarseAggregate    0.018086   0.009392   1.926 0.054425 .  
## FineAggregate      0.020190   0.010702   1.887 0.059491 .  
## Age                0.114222   0.005427  21.046  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.4 on 1021 degrees of freedom
## Multiple R-squared:  0.6155, Adjusted R-squared:  0.6125 
## F-statistic: 204.3 on 8 and 1021 DF,  p-value: < 2.2e-16

Посмотрим на корреляцию:

corrplot(cor(concr.data))

xyplot(concr.data$Water ~ concr.data$Superplasticizer, panel=panel)

Построим модель, откинув не сильно значимые признаки агрегаторов, а также добавим Water:Superplasticizer, т.к. между ними явно есть зависимость

model1 <- lm(ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag + FlyAsh + Superplasticizer + Age + Water:Superplasticizer, data=concr.data)
summary(model1)

## 
## Call:
## lm(formula = ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag + 
##     FlyAsh + Superplasticizer + Age + Water:Superplasticizer, 
##     data = concr.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.439  -6.400   0.209   6.940  33.095 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            45.502486   5.121626   8.884  < 2e-16 ***
## Cement                  0.104717   0.004189  25.001  < 2e-16 ***
## Water                  -0.304965   0.026125 -11.673  < 2e-16 ***
## BlastFurnaceSlag        0.082116   0.004967  16.532  < 2e-16 ***
## FlyAsh                  0.048647   0.008447   5.759 1.12e-08 ***
## Superplasticizer       -2.351854   0.477719  -4.923 9.93e-07 ***
## Age                     0.120990   0.005502  21.989  < 2e-16 ***
## Water:Superplasticizer  0.015927   0.002890   5.511 4.52e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.26 on 1022 degrees of freedom
## Multiple R-squared:  0.6252, Adjusted R-squared:  0.6226 
## F-statistic: 243.5 on 7 and 1022 DF,  p-value: < 2.2e-16

Немного лучше, но еще не очень.

Посмотрим на графики распределений предикторов:

marginal.plot(concr.data)

Age выглядит совсем плохо, попробуем прологарифмировать

marginal.plot(log(concr.data$Age))

xyplot(concr.data$ConcreteCompressiveStrength ~ log(concr.data$Age), panel=panel)

Стало лучше Посмотрим модель

model2 <- lm(ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag + FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer, data=concr.data)
summary(model2)

## 
## Call:
## lm(formula = ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag + 
##     FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer, 
##     data = concr.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.9591  -4.4094   0.0202   4.2865  28.9493 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            19.464258   3.468202   5.612 2.57e-08 ***
## Cement                  0.109381   0.002929  37.343  < 2e-16 ***
## Water                  -0.283281   0.017676 -16.027  < 2e-16 ***
## BlastFurnaceSlag        0.083428   0.003473  24.023  < 2e-16 ***
## FlyAsh                  0.050869   0.005905   8.614  < 2e-16 ***
## Superplasticizer       -1.164789   0.325079  -3.583 0.000356 ***
## log(Age)                8.712859   0.192143  45.346  < 2e-16 ***
## Water:Superplasticizer  0.007442   0.001965   3.787 0.000161 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.177 on 1022 degrees of freedom
## Multiple R-squared:  0.8167, Adjusted R-squared:  0.8154 
## F-statistic: 650.4 on 7 and 1022 DF,  p-value: < 2.2e-16

Намного лучше, но Мне все еще не нравится изломанность графика воды Поэтому можно попробовать добавить фактор для воды:

concr.data$Water.factor <- concr.data$Water < 225
model3 <-lm(ConcreteCompressiveStrength ~ Cement + Water.factor:Water + BlastFurnaceSlag + FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer, data=concr.data)
summary(model3)

## 
## Call:
## lm(formula = ConcreteCompressiveStrength ~ Cement + Water.factor:Water + 
##     BlastFurnaceSlag + FlyAsh + Superplasticizer + log(Age) + 
##     Water:Superplasticizer, data = concr.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.7822  -4.3935   0.0625   4.5636  26.9560 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             36.895894   4.578324   8.059 2.14e-15 ***
## Cement                   0.106611   0.002925  36.446  < 2e-16 ***
## BlastFurnaceSlag         0.080584   0.003456  23.316  < 2e-16 ***
## FlyAsh                   0.042757   0.005986   7.143 1.74e-12 ***
## Superplasticizer        -2.207334   0.368434  -5.991 2.89e-09 ***
## log(Age)                 8.471278   0.193891  43.691  < 2e-16 ***
## Water.factorFALSE:Water -0.334456   0.019574 -17.087  < 2e-16 ***
## Water.factorTRUE:Water  -0.368904   0.022963 -16.065  < 2e-16 ***
## Water:Superplasticizer   0.013515   0.002207   6.122 1.31e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.069 on 1021 degrees of freedom
## Multiple R-squared:  0.8224, Adjusted R-squared:  0.821 
## F-statistic: 590.8 on 8 and 1021 DF,  p-value: < 2.2e-16

Немного поперебирав значение порога получаем высоко значимый признак Water.factorTrue:Water, который говорит как влияет количство воды до 225 на прочность

Попробуем также и возраст сделать фактором:

concr.data$Age.factor <- concr.data$Age < 150
model4 <-lm(ConcreteCompressiveStrength ~ Cement + Water.factor:Water + BlastFurnaceSlag + FlyAsh + Superplasticizer + Age.factor:log(Age) + Water:Superplasticizer, data=concr.data)
summary(model4)

## 
## Call:
## lm(formula = ConcreteCompressiveStrength ~ Cement + Water.factor:Water + 
##     BlastFurnaceSlag + FlyAsh + Superplasticizer + Age.factor:log(Age) + 
##     Water:Superplasticizer, data = concr.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.7246  -4.3742  -0.1258   4.3747  26.5782 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              32.336381   4.433372   7.294 6.04e-13 ***
## Cement                    0.107538   0.002816  38.185  < 2e-16 ***
## BlastFurnaceSlag          0.078804   0.003331  23.657  < 2e-16 ***
## FlyAsh                    0.040766   0.005764   7.073 2.82e-12 ***
## Superplasticizer         -1.907606   0.356011  -5.358 1.04e-07 ***
## Water.factorFALSE:Water  -0.306365   0.019084 -16.054  < 2e-16 ***
## Water.factorTRUE:Water   -0.353940   0.022154 -15.976  < 2e-16 ***
## Age.factorFALSE:log(Age)  7.331764   0.224633  32.639  < 2e-16 ***
## Age.factorTRUE:log(Age)   9.301197   0.207623  44.799  < 2e-16 ***
## Water:Superplasticizer    0.011394   0.002137   5.333 1.19e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.801 on 1020 degrees of freedom
## Multiple R-squared:  0.8357, Adjusted R-squared:  0.8343 
## F-statistic: 576.5 on 9 and 1020 DF,  p-value: < 2.2e-16

И это тоже улучшило модель

Проверим модели

anova(model1, model2, model3, model4)

## Analysis of Variance Table
## 
## Model 1: ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag + 
##     FlyAsh + Superplasticizer + Age + Water:Superplasticizer
## Model 2: ConcreteCompressiveStrength ~ Cement + Water + BlastFurnaceSlag + 
##     FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer
## Model 3: ConcreteCompressiveStrength ~ Cement + Water.factor:Water + BlastFurnaceSlag + 
##     FlyAsh + Superplasticizer + log(Age) + Water:Superplasticizer
## Model 4: ConcreteCompressiveStrength ~ Cement + Water.factor:Water + BlastFurnaceSlag + 
##     FlyAsh + Superplasticizer + Age.factor:log(Age) + Water:Superplasticizer
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   1022 107645                                  
## 2   1022  52647  0     54998                     
## 3   1021  51013  1      1633 35.316 3.843e-09 ***
## 4   1020  47178  1      3835 82.924 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Последняя модель самая лучшая

xyplot(residuals(model1) ~ fitted(model1), panel=panel)

xyplot(residuals(model2) ~ fitted(model2), panel=panel)

xyplot(residuals(model3) ~ fitted(model3), panel=panel)

xyplot(residuals(model4) ~ fitted(model4), panel=panel)

Более менее линейно

cv <- function(model) { 
  tune(lm, model$call$formula, data = concr.data, tunecontrol = tune.control(sampling = "cross")) 
}
cv(model1)

## 
## Error estimation of 'lm' using 10-fold cross validation: 107.0302

cv(model2)

## 
## Error estimation of 'lm' using 10-fold cross validation: 52.16315

cv(model3)

## 
## Error estimation of 'lm' using 10-fold cross validation: 50.36774

cv(model4)

## 
## Error estimation of 'lm' using 10-fold cross validation: 46.80454

На убывание, однако итоговая ошибка все-равно очень велика относительно диапозона прочности